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Simulating individual-based models of bacterial chemotaxis with 
asymptotic variance reduction 



Mathias Rousset* and Giovanni Samaey^ 

We discuss variance reduced simulations for an individual-based model of chemotaxis of 
bacteria with internal dynamics. The variance reduction is achieved via a coupling of this 
model with a simpler process in which the internal dynamics has been replaced by a direct 
gradient sensing of the chemoattractants concentrations. In the companion paper , we 
have rigorously shown, using a pathwise probabilistic technique, that both processes 
converge towards the same advection-diffusion process in the diffusive asymptotics. In 
this work, a direct coupling is achieved between paths of individual bacteria simulated 
by both models, by using the same sets of random numbers in both simulations. This 
coupling is used to construct a hybrid scheme with reduced variance. We first compute 
a deterministic solution of the kinetic density description of the direct gradient sensing 
model; the deviations due to the presence of internal dynamics are then evaluated via 
the coupled individual-based simulations. We show that the resulting variance reduction 
is asymptotic, in the sense that, in the diffusive asymptotics, the difference between the 
two processes has a variance which vanishes according to the small parameter. 

MSG: 35Q80, 92B05, 65C35. 
Key Vifords: asymptotic variance reduction, bacterial chemotaxis, velocity-jump 
process. 

1. Introduction 

The motion of flagellated bacteria can be modeled as a sequence of run phases, 
during which a bacterium moves in a straight line at constant speed. Between 
two run phases, the bacterium changes direction in a tumble phase, which takes 
much less time than the run phase and acts as a reorientation. To bias movement 
towards regions with high concentration of chemoattractant, the bacterium adjusts 
its turning rate by increasing, resp. decreasing, the probability of tumbling when 
moving in an unfavorable, resp. favorable, direction Since many species are 
unable to sense chemoattractant gradients reliably due to their small size, this 
adjustment is often done via an intracellular mechanism that allows the bacterium 
to retain information on the history of the chemoattractant concentrations along its 
path's!. The resulting model, which will be called the "internal state" or "fine-scale" 
model in this text, can be formulated as a velocity-jump process, combined with 
an ordinary differential equation (ODE) that describes the evolution of an internal 
state that incorporates this memory effect 
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2 Variance reduced simulation of chemotaxis 

The probability density distribution of the velocity-jump process evolves accord- 
ing to a kinetic equation, in which the internal variables appear as additional dimen- 
sions. A direct deterministic simulation of this equation is therefore prohibitively 
expensive, and one needs to resort to a stochastic particle method. Unfortunately, 
a direct fine-scale simulation using stochastic particles presents a large statistical 
variance, even in the diffusive asymptotic regime, where the behavior of the bacte- 
rial density is known explicitly to satisfy a Keller-Segel advection-diffusion equation. 
Consequently, it is difficult to assess accurately how the solutions of the fine-scale 
model differ from their advection-diffusion limit in intermediate regimes. We refer 
toUmZIIIlfor numerous historical references on the Keller-Segel equation, toEElfor 
formal derivations (based on moment closures) of convergence of the the velocity- 
jump process with internal state to the Keller-Segel equation, and to the companion 
paper 1211 for a proof of this convergence based on probabilistic arguments. 

In this paper, we propose and analyze a numerical method to simulate 
individual-based models for chemotaxis of bacteria with internal dynamics with 
reduced variance. The variance reduction is based on a coupling technique (con- 
trol variate): the main idea is to simultaneously simulate, using the same random 
numbers, a simpler, "coarse" process where the internal dynamics is replaced by 
a direct "gradient sensing" mechanism (see Pp3|25 | j^j. references on such gradient 
sensing models). The probability density of the latter satisfies a kinetic equation 
without the additional dimensions of the internal state, and converges to a similar 
advection-diffusion limit, see e.g. |6|i5|24|26 [ 

More precisely, we consider the coupling of two velocity-jump processes with 
exactly the same advection-diffusion limit (for both processes, a proof of the latter 
convergence has been obtained in the companion paper . The first one is the 



process with internal dynamics described by the system of equations (2.4 1, with 



jump rates satisfying (2.3 1 and internal dynamics satisfying the two main assump- 



tions (2.7) and (2.9) as well as other technical assumptions contained in Section 2.2 



The second one is the process with direct gradient sensing (2.12), with jump rates 



satisfying (2.14) with (2.15). The same random numbers are used for the two pro- 
cesses in the definition of both the jump times and the velocity directions. The main 
contributions of the present paper are then twofold : 



From a numerical point of view, we couple two systems of many parti- 
cles consisting of different realizations of the fine-scale (internal state) and 
coarse (gradient sensing) processes, simulated with similar discretization 



schemes (see Section 3.1) and the same random numbers. Additionally, for 
the coarse model, also the continuum description is simulated on a spatial 
grid. Evolution of the fine-scale model is then evaluated on the grid at hand 
by adding the difference between the two coupled particle descriptions to 
the evolution of the coarse continuum description. This can be seen as us- 
ing the coarse process as a control variate. As coupling is lost over time, 
a regular reinitialization of the two coupled particle systems is included in 
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Model 


Numerics 


Internal state model-Individual description 
Gradient sensing model-Individual description 
Gradient sensing model-Density description 


Monte-Carlo method 
Monte-Carlo method (coupled) 
Grid method 



Table 1. Models and numerical methods used in the variance reduction algorithm. 



the algorithm. The asymptotic variance reduction method can be depicted 
usmg Table [l]( see also Figure [l] in Section |3| . 

From an analysis point of view, we rigorously prove, using probabilistic 
arguments, L^-bounds on the position difference at diffusive times of both 
processes in terms of appropriate powers of the scaling parameter (The- 
orem 4.7). This analysis shows that the numerical variance reduction ob- 



tained by the method is asymptotic (the remaining statistical variance van- 
ishes with the small scaling parameter), and provides upper bounds on the 
rate of the latter convergence. We argue in Section |4.3| that this bound 
should be sharp in relevant regimes. 



The idea of asymptotic variance reduction is a general and very recent idea 
in scientific computing that appears when using hybrid Monte-Carlo/PDE (partial 
differential equation) methods. While the use of control variates is fairly common in 
Monte Carlo simulation, the only explicit attempt to develop asymptotic variance 
reduction in hybrid methods, up to our knowledge, can be found in and related 
papers in the context of the Boltzmann equation, and is based on an importance 
sampling approach (see also and related papers). 

This paper is organized as follows. In Section [2j we briefly review the mathe- 
matical models that we will consider, and summarize the results on their advection- 
diffusion limit that were obtained in the companion paper t^H. In Section [sj we 
proceed to describe the coupled numerical method in detail. Section |4] is dedicated 
to a detailed analysis of the variance reduction of the coupled method. The main 
result is the following: up to a time shift perturbation, the variance of the difference 
between the two coupled process scales according to the small parameter of the 
diffusive asymptotics. We proceed with numerical illustrations in Sections [5] and [6] 
and conclude in Section [7] with an outlook to future research. 



2. Particle-based models for bacterial chemotaxis 

In this paper, we directly present the models of interest in a nondimensional form, 
incorporating the appropriate space and time scales. For a discussion on the chosen 
scaling, we refer to the companion paper I^Sl. 
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4 Variance reduced simulation of chemotaxis 

2.1. Model with internal state 

We consider bacteria that are sensitive to the concentration of m chemoattractants 
{Pi{x))^i, with Pi{x) >0 for xgM''. While we do not consider time dependence of 
chemoattractant via production or consumption by the bacteria, a generahzation 
to this situation is straightforward, at least for the definition of the models and 
the numerical method. Bacteria move with a constant speed v (run), and change 
direction at random instances in time (tumble), in an attempt to move towards 
regions with high chemoattractant concentrations. As in we describe this be- 
havior by a velocity-jump process driven by some internal state y G Y C M" of each 
individual bacterium. The internal state models the memory of the bacterium and 
is subject to an evolution mechanism attracted by a function -0 : M'" M" of the 
chemoattractants concentrations, 

where x is the present position of the bacterium. S is assumed to be smooth with 
bounded derivatives up to order 2. We refer to Section 2.2 in ^ for some typical 
choices; in the concrete example in this paper, we choose m = n = l, and S{x) = 
Pi{x). By convention, the gradient of S{x) is a matrix with dimension 

V5(x)eM"^''. 

We then denote the evolution of each individual bacterium position by i i— > X(, with 
normalized velocity 

at 

with S'^~^ the unit sphere in R''. Hence, Vt represents the direction and the pa- 
rameter e represents the size of the velocity. The evolution of the internal state is 
denoted by ti—>-Yt. The internal state adapts to the local chemoattractant concen- 
tration through an ordinary differential equation (ODE), 

-^=F,{YuS{X,)), (2.1) 

which is required to have a unique fixed point y* ^ S{x*) for every fixed value x* G 
M.'^. We also introduce the deviations from equilibrium z = S{x) — y. The evolution 
of these deviations is denoted as 

t^Zt = S{Xt)~Yt. 

The velocity of each bacterium is switched at random jump times ('rn)ri>i that are 
generated via a Poisson process with a time dependent rate given by \{Zt), where 
z^\{z) is a smooth function satisfying 

0<A„,in<A(z)<A„iax, (2.2) 

as well as, 

X{z) = \o-b^z + cxO(\zt), (2.3) 
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with h e M", /c > 2, and ca > is used to keep track of the non-linearity in the anal- 
ysis. The new velocity at time T„ is generated at random according to a centered 
probability distribution M.{dv) with j vM.{dv) — Q^ typically 

Ai{dv) = cTgd-i (dw), 

where crgd-i is the uniform distribution on the unit sphere. 

The resulting fine-scale stochastic evolution of a bacterium is then described by 
a left continuous with right limits (Icrl) process 

which satisfies the following differential velocity-jump equation : 



At 

m 



= F,{YuS{Xt)), 



(2.4) 



w\t\iZt:^S{Xt)-Yt, 



Vt^Vn foriG[r„,r„+i] , 

with initial condition Xo,Vo eM"*, Fq eK" and Tq = 0. In ([O]), {9n)n>i denote i.i.d. 
random variables with normalized exponential distribution, and (V„)^>-^ denote 
i.i.d. random variables with distribution Ai{dv). 

Example 2.1. For concreteness, we provide a specific example, adapted fromlS]^ 
which will also be used later on to illustrate our results numerically. We consider 
m= 1, i.e., there is only one chemoattractant S{x) and the internal dynamics (2.1 1 
reduces to a scalar equation 



<^y[t) _ S{x)-y _ z 

At T T 



(2.5) 



For the turning rate A(z), we choose the following nonlinear strictly decreasing 
smooth function 



A(z) = 2Ao 





1 / 






— arctan | 


v2Ao^^)) 




TT ' 





(2.6) 



2.2. Asymptotic regimes and linearizations 



In the adimensional models ( 2.4 )-( 2.12), a small parameter e is introduced, which 
can be interpreted as the ratio of the typical time of random change of velocity 
direction to the typical time associated with chemoattractant variations as seen by 
the bacteria. In the present section, we present the main scaling assumptions on 
the internal dynamics that are necessary to consider the asymptotic regime e^l, 
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as well as the corresponding notation. Most of these assumptions were introduced 
in the companion paper to which we refer for more details and motivation. 

Remark 2.1 (Notation). Throughout the text, the Landau symbol O denotes a 
deterministic and globally Lipschitz function satisfying 0(0) = 0. Its precise value 
may vary from line to line, may depend on all parameters of the model, but its 
Lipschitz constant is uniform in e. In the same way, we will denote generically by 
C > a deterministic constant that may depend on all the parameters of the model 
except e. 



We assume that the ODE (2.1) driving the internal state is well approximated 



by a near equilibrium evolution equation in following sense : 
Assumption 2.2. We have 

F,{y,s) = -T-\y-s) + e'-'cFO{\s-y\^), 
where S>0 and t, e M"^" is an invertible constant matrix. 



(2.7) 



The constant cp is used to keep track of the dependance on the non-linearity of 
in the analysis. For technical reasons that are specific to the analysis of the time 



shift in the coupling (Lemma 4.6), we need the following assumption on the scale of 
the non-linearity (this assumption was not necessary in the companion paper : 

Assumption 2.3. Assume that jy — s| —0{€^), then there is a 7>0 such that: 

\{y-s)-T,F,{y,s)\^CFO{e'+''). (2.8) 

We now formulate the necessary assumptions on the scale of the matrix : 

Assumption 2.4. There is a constant C> such that for any t>0, one has 

<Ce-*''"'/^. (2.9) 

Assumption 2.5. There is a constant C > such that 



sup 



<Ce-\ 



Assumption |2.4| is used to ensure exponential convergence of the linear ODE with 
time scale at least of order e^"^ . Assumption |2.5| is necessary for technical reasons 

(which is given in 



IS symmetric 



in the proof of Lemma 4.5 (which is given inl^. Note that if r^^^ 

and strictly positive with a lower bound of order 0{€^~^) on the spectrum, then 

Assumption [2^ and Assumption [23] are satisfied. This implies that the slowest time 

s 



scale of the ODE (2.1) is at worse of order 

-1 



Assumption 2.4 implies ||r< 



Indeed, integration over time in 
Moreover, we assume that the solution of the 



ODE driving the internal state in (2.4 1 satisfies the following long time behavior 
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Assumption 2.6. Consider a path tt-^St such that suptg^^^ 

assume \Yo — So\=0{e^). Then the solution of 

dYt 



dSt 



dt 



:0(e), and 



dt 



satisfies 



sup \Yt~St\ 
te[o,+oo] 



Assumption |2.6| is motivated by the case of hnear ODEs; indeed, in that case. 



Assumption 2.6 is a consequence of Assumption 2.4 (see Section 3.3 inlSSI). Recalhng 
that Zt = St — Yt, we will assume in the remainder of the paper that the solution 
of (2.4) satisfies. 



\Zo\=0{e'), 



as well as 



Then, under Assumption |2.2[ Assumption |2.4[ Assumption |2.5[ and Assumption |2.6| 
we obtain the bounds (see Section 3.3 in 



sup 

ie[0,+oo] 



r.-^e- 



\St-Yt) =0(1), 



(2.10) 



sup {\Zt\)=0{e'), 

*G[0,+oo] 



as well as 



sup 

tG[0,+O( 



'Z, 



= 0{l). 



Finally, we will make use of the following notation. Denoting eg = || V^S*]!^ , 
with the Hessian, the error terms due to non-linearities of the problem will be 
handled through 

fnl(e) ■.^CFe^ + cse + c^e''^-^=o{\), 
1nl2(e):=Cf(e* + e^)+C5(e + e*')+CAe*^-'^-i=o(l). 



(2.11) 



Note that nl(e)<nl2(e) and nl(e) --nl2(e) \i 5<l and 5<^. RecaU that fc>2 is 
defined in p^] ). 

Note that all the assumptions of the present section hold in the following case: (i) 
the ODE is linear (i^^ is linear), (ii) the associated matrix is symmetric positive 
and satisfies 

for some (5 > 0, (iii) the initial internal state is close to equilibrium in the sense that 
Te^^Zo = 0(1). The assumptions of the present section can be seen as technical 
generalizations to non-linear ODEs with non-symmetric linear part. 
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2.3. Model with direct gradient sensing (control process) 



We now turn to a simplified, "coarse" model, in which the internal process (2.1 ), and 



the corresponding state variables, are eliminated. Instead, the turning rate depends 
directly on the chemoattractant gradient. This process will be called the control 
process since it will be used in Section [3] as a control variate to perform variance 
reduced simulations of (12. 41). 



The control process is a Markov process in position-velocity variables 
which evolves according to the following differential velocity-jump equations : 



dt 



KiX'^^V^dt- 



(2.12) 



with initial condition Xq, Vq G. 



= v„ forteK,r„%i] , 

¥K In (|2.12|), 



n>l 



denote i.i.d. random variables 
with normalized exponential distribution, and {yn)n>i denote i.i.d. random vari- 
ables with distribution Ai{dv). The turning rate of the control process is assumed 
to satisfy 

0<An,i„<AJ(x,t;)<A,„ax, (2.13) 

as well as 

K{x,v):=\a-eAj{x)v + 0{e^). 



(2.14) 



Typically, \'^{x,v) is a function of V5(x), so that the model (2.14) may describe 



a large bacterium that is able to directly sense chemoattractant gradients. When 
m = n = l, and the turning rate (2.14) is proportional to VS'(a;)uGM, it can be in- 



terpreted as follows: the rate at which a bacterium will change its velocity direction 
depends on the alignment of the velocity with the gradient of the chemoattractant 
concentration VS'(a;), resulting in a transport towards areas with higher chemoat- 
tractant concentrations. In the companion paper it is shown that, when the 
vector field {x) £ M'* is given by the formula 

M^)^b^TT^^Six), (2.15) 
asymptotic consistency with the internal state model is obtained in the limit when 



e— ^0 (see Section 2.4 1 



2.4. Kinetic formulation and advection- diffusion limits 

We now turn to the kinetic description of the probability distribution of the pro- 
cesses introduced above. The probability distribution density of the fine-scale pro- 
cess with internal state at time t with respect to the measure dxM.{dv)dy is de- 
noted as p{x,v,y,t), suppressing the dependence on e for notational convenience, 
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and evolves according to the Kolomogorov forward evolution equation (or master 
equation). In the present context, the latter is the following kinetic equation 



dtp + ev ■ V^p + divy (F, ix,y)p) ^ X {S {x) - y) {R{p) - p) , 



(2.16) 



where 



i?(p):=/ pi;v,-)Midv) 

is the operator integrating velocities with respect to Ai, and F{x,y) and X{z) are 
defined as in Section |2.1| Similarly, the distribution density of the control process 
is denoted as p'^(x,w,t), and evolves according to the kinetic/master equation: 

dtp' + ev ■ V,p^ = (i?(AJp) - AJp) , (2.17) 

with defined as in (2.14). The reader is referred to for the derivation of master 



equations associated to Markov jump processes. 

In the companion paper we show, using probabilistic arguments, that, in 



the limit e— 7>0, both the equation for the control process (2.17) and the equation 



for the process with internal state (2.16) converge to an advection-diffusion limit 



on diffusive time scales. Convergence is to be understood pathwise, in the sense of 
convergence of probability distribution on paths endowed with the uniform topology. 
Denote diffusive times by 

and the processes considered on diffusive time scales as 

For the control process, we then have the following theorem : 



Theorem 2.1. For e— )-0, the process ti-^X^' , solution of (2.12), converges to- 
wards an advection-diffusion process, satisfying the stochastic differential equation 
(SDE) 



Ao 



di- 



2D 



1/2 



dW, 



(2.18) 



where 1 1-> is a standard Brownian motion, the parameters Xq and Aq := limg^yo ^£ 



originate from the turning rate (2.14), and the diffusion matrix is given by the 



covariance of the Maxwellian distribution : 



D- 



v®vM{dv)&M!^'"^. 



(2.19) 



In particular, this result implies that, at the level of the Kolomogorov/master 
evolution equation, the evolution of the position bacterial density. 



n'''\x,t):=n''{x,t/e^):^ / p''{x,v,t/e^)M{dv) 



(2.20) 
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converges to the advection-difFusion equation 

1 



-^0 



(2.21) 



on diffusive time scales as e— >0. 

In the same way, a standard probabilistic diffusion approximation argument can 



be used to derive the pathwise diffusive limit of the process with internal state ( 2.4 ) : 



Theorem 2.2. For e— >0, the process ti-^X^, solution of (2.4), converges to 



wards an advection-diffusion process, satisfying the stochastic differential equation 

-VS{x), 



(SDE) (2.18), where Aq originates from 

Aa{x)^h^ lim 



£^0 AoTf +Id 



(2.22) 



in which, b, Tj, and Aq were introduced in (2.7)-( |2.3| as parameters of the process 
with internal state, and IdeM"^" is the identity matrix. Again, the diffusion matrix 



D is given by the covariance of the Maxwellian distribution (2.19) 



Also here, convergence needs to be understood in terms of probability distribu- 
tion on paths endowed with the uniform topology. 

Introducing the bacterial density of the process with internal state as 



n{x,t) : 



p{x,v,y,t)M{dv)dy, 



(2.23) 



this implies that the evolution of n converges to (2.21) on diffusive time scales in 
the limit of 0. 



3. Numerical method 



To simulate the process with internal state, solving the kinetic equation (2.16) 



over diffusive time scales can be cumbersome, due to the additional dimensions as- 
sociated with the internal state. The alternative is to use to stochastic particles. 



However, a particle-based simulation of equation (2.161 is subject to a large statis- 
tical variance of the order ©(TV"^/^), where N is the number of simulated particles. 
The asymptotic analysis shows that the position bacterial density approaches an 
advection-diffusion limit (2.21) when e^O. Consequently, to accurately assess the 



deviations of the process with internal state (2.4) as compared to its advection- 



diffusion limit (for small but non- vanishing (intermediate) values of e), the required 
number of particles needs to increase substantially with decreasing e, which may 
become prohibitive from a computational point of view. 

In this section, we therefore propose a hybrid method, based on the principle of 
control variates, that couples the process with internal dynamics with the control 
process, which is simulated simultaneously using a grid-based method. We first 
describe the variance reduction technique (Section 3.1). The analysis in Section |4] 
will reveal that the variance reduction is asymptotic, in the sense that the variance 
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vanishes in the diffusion limit. To ensure this asymptotic variance reduction during 
actual simulations, one needs to ensure that the time discretization preserves the 
diffusion limits of the time-continuous process. An appropriate time discretization 
is discussed in Section [3^ 

3.1. Coupling and asymptotic variance reduction 

The proposed variance reduction technique is based on the introduction of a control 
variate that exploits a coupling between the process with internal state and the con- 
trol process. While the idea of control variates for Monte Carlo simulation is already 
well known, see e.g. and references therein, the coupling that is proposed here is 
particular (we call it asymptotic), since the difference of the coupled processes, and 
hence the variance, vanishes with an estimable rate in the diffusion limit e— >0, as 
will be shown in Section |4| 



3.1.1. The control variates 

Let us first assume that we are able to compute the exact solution of the kinetic 



equation for the control process, (2.171, with infinite precision in space and time. 

The algorithm of asymptotic variance reduction is based on a coupling between 
an ensemble of realizations evolving according to the process with internal state 



(2.4), denoted as 

and an ensemble of realization of the control process (2.121, denoted as 

We denote the empirical measure of the particles with internal state in position- 
velocity space as 

1 ^ 

i—l 

and, correspondingly, the empirical measure of the control particles as 

C.N _ 1 r 
i—l 

A coupling between the two ensembles is obtained by ensuring that both sim- 
ulations use the same random numbers {9„)n>i and (V„)„>o, which results in a 
strong correlation between {XI, V^) and (X^''^, Vj'^'*) for each realization. Simulta- 
neously, the kinetic equation for the control process ( 2.17| ) is also solved using a 



deterministic method (which, for now, is assumed to be exact). We formally denote 
the corresponding semi-group evolution as 

e*^', with L"(p") = -ewV:,p" + (i?(A^p'=)-AJp'=). 
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Besides the two particle measures and A*^'^, we denote by JI^ the variance 
reduced measure, which will be defined by the algorithm below. Since, with increas- 
ing diffusive time, the variance of the algorithm increases due to a loss of coupling 
between the particles with internal state and the control particles, the variance 
reduced algorithm will also make use of a reinitialization time step Stj-i, which is 
defined on the diffusive time scale. The corresponding time instances are denoted 
as tn = nStri on the diffusive time scale, or equivalently, on the original time scale 
as tn = nStri/e^. 

Starting from an initial probability measure /io at time t~0, we sample fiQ to 
obtain the ensemble {Xl,Vf ,Yf^ corresponding to fiQ , and then set := , 
i.e., Xq'^ = Xq and Vq''' — Vq for all i = l,...,N. Furthermore, we set the variance 
reduced estimator as JLq := /zq =E(/i^). We then use the following algorithm to 
advance from i„ to (see also Figure (fTl) 



Algorithm 3.1. At time we have that the particle measure /i^' — fJ-^ , and the 
variance reduced measure is given by Jl^ . To advance from time i„ to we 
perform the following steps : 



Evolve the particles {Xl,Vf j from <„ to tn+i, according to (2.4), 



Evolve the particles lATj''^, V/'*^! according to (2.121, using the same 



random numbers as for the process with internal state, 

• Compute the variance reduced evolution 

• Reinitialize the control particles by setting 

k:+i-xi^i, z=i,...,7v, 

i.e., we set the state of the control particles to be identical to the state of 
the particles with internal state. 



In formula (3.1), we use the symbol t~_^-^ to emphasize that the involved parti- 
cle positions and velocities are those obtained before the reinitialization. An easy 
computation shows that the algorithm is unbiased in the sense that for any n>0, 

E(/2fJ=E(/xfJ, 

since the particles with internal dynamics are unaffected by the reinitialization, and, 
additionally, 

E/x5;^^=EM^e^"/^'^^ 

Moreover, the variance is controlled by the coupling between the two processes. 
Indeed, using the independence of the random numbers between two steps of Algo- 



rithm 3.1 and introducing ip as a position and velocity dependent test function) 
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Variance reduction with reinitiaiization 



- Particies with internai variabies 
Particies with gradient sensing 

Deterministic grid-based simuiation: gradient sensing 

- Variance reduced simulation 




Fig. 1. A schematic description of Algorithm |3.1| The dashed Hne represent the evolution of 
N bacteria with internal state. The dotted line represent the coupled evolution of A'^ bacteria 
with gradient sensing, subject to regular reinitializations. The dashed-dotted line is computed 
according to a deterministic method simulating the density of the model with gradient sensing, and 
subject to regular reinitializations. The solid line is the variance reduced simulation of the internal 
state dynamics, and is computed by adding the difference between the particle computation with 
internal state, and the particle simulation with gradient sensing to the deterministic gradient 
sensing simulation. At each reinitialization step, the two simulations (deterministic and particles) 
of the gradient sensing dynamics are reinitialized to the values of their internal state simulation 
counterpart (as represented by the arrows). 



we get (according to the main result of Theorem |4.7| 



var(^f^ M) = J2^L^{^) ~ 



fc=i 
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and thus 



var(7I-M)<Cn^^^^±^, (3.2) 



where in the last hne, C is independant of n, e, and N. 



In some generic situations (see Section 4.3), we can argue that the statistical 
error in Algorithm |3.1| coming from the coupling is "sharp" with respect to the 
order in e. This means that the difference between the probabihty distribution of 
the model with internal state and the probability distribution of the model with 
gradient sensing is of the same order. This would imply that, with the asymptotic 
variance reduction technique, one is able to reliably assess the true deviation of the 
process with internal variables from the control process using a number of particles 
N that is independent of e. 



3.1.2. Time-space discretization of the kinetic equation 

For both processes, the probability density distribution of the position and velocity 
can then be computed via binning in a histogram, or via standard kernel density 
estimation EZEH using a kernel Khi in which the chosen bandwidth h can be based 
on the grid used for the deterministic simulation and on the data, for instance using 
the Silverman heuristic f^Sl. This yields 

1 ^ 

p^(x,«,f) = — ^i^;,(a:-X^«-F/), (3.3) 
i=l 

and similarly 

1 ^ 

p'k{x,v,t)^—Y,K,{x~Xr,v-Vn. (3.4) 

The important point is that the solution p'^{x,v,t) for the control process may 
be accurately approximated with a deterministic grid-based method, which ensures 
that 

\p^x,v,i/e')-pl,,^ix,v,t/e^)\^0{Sx') + OiSt'^), 
for some integers k,l>l. 

Remark 3.1. In this text, we will perform the simulations in one space dimension 
using a third-order upwind-biased scheme, and perform time integration using the 
standard fourth order Runge-Kutta method, i.e., 1 = 3 and fc = 4. 

Then, the unbiased nature of the variance-reduced estimator is conserved up to 
0{Sx'')+0{St'') + 0{h) discretization errors. 
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3.2. Time discretization of velocity-jump processes 



When simulating equation (2.4), a time discretization error originates from the fact 
that the equation for the evolution of the internal state, and hence the evolution 
of the fluctuations Zt, is discretized in time. This results in an approximation of 
the jump times (T„)„>i, and hence of Xt- To retain the diffusion limits of the 
time-continuous process, special care is needed. We now briefly recall the time 
discretization procedure that was proposed and analyzed in the companion paperl^. 



For ease of exposition, we consider the scalar equation (2.51 for the internal state; 
generalization to nonlinear systems of equations is briefly discussed in 



Consider flrst the linear turning rate (2.3 1. In that case, we define a numerical 



solution (Xf as follows. Between jumps, we discretize the simulation in 

steps of size St and denote by (Xf^^^^Zf^*,^) the solution at tn,k=T^* + k5t. The 
numerical solution for \tn,k,tn,k+i] is given by 

fxf =X5 + eV„(i-<„,fe) 

\zf =exp(-(t-t„.fc)r,-i)Z^*^+er,(ld-exp(-(t-t„,fc)T,-i))V5(X;^\)V„. 

(3.5) 

We denote by K>Q the integer such that the simulated jump time T^\i^ 

[tn,K^'tn,K+i\- To find T^\n we first approximate the integral /^^+^ A(Zf*)di us- 
ing 



/ A(Zf )dt= ^ / "''^'A(Zf )dt-H / A(Zf )dt, (3.6) 
and then compute : 

^ A(Zf )di = Ao(5i-&^ (id-e-**^'"') r^Z^^ 



eb' [5tT,-{lA-e~''^'~ W]^S{Xi\)Vn. (3.7) 



The jump time T^f^^^ can then be computed as the solution of 

A(zf )di = 0„+i - E / " Hz!')dt, (3.8) 
using a Newton procedure. It is shown in that the results on the diffusion limit 



as outlined in Section 2.4 1 are not affected by the discretization. 



We now consider a general nonlinear turning rate. We again discretize in time 

t 

-1-1 



to obtain the time-discrete solution (3.5 1. The jump time T^+i is now computed by 



linearizing (2.6) in each time step. 



A** (Zf , = MZf^,,) + (Zf * - <,) . (3.9) 
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One then replaces equation (3.6 1 by 



and proceeds in the same way as for the hnear case. Also in this case, the diffusive 
limit is recovered in an exact fashion for the time-discretized process. 

Remark 3.2 (Discretization of the control process). For the control process, 
there is no internal state. Hence, the only time dependence of turning rate is due 
to the spatial variation of V S{x), which can be treated by discretizing the integral 
of A in the same way as above. 

4. Asymptotic variance reduction of the coupling 

In this section, we show how the difference between the two coupled processes on 
diffusive time scales t — t/e^ behaves in the limit of e— >0. We first recall some 



notation from the companion paper (Section 4.1 1, after which we state and prove 
the main theorem (Section |4.2|. 



4.1. Notations and asymptotic estimates 



The main theorem in Section [4.2| relies on following auxiliary definitions and lemmas 
that were given in the companion paper 

Definition 4.1. We denote by m:M^-M"^" the function 

m(<):=ir,- (ld-e-*^="')r,2, (4.1) 



whose derivative is given by 

m'(t)=T, (id -e-*^= 
This function satisfies the following lemmas : 
Lemma 4.2. Let 9 be an exponential random variable of mean 1. Then: 



see equation (2.15) 



Lemma 4.3. For allt^R, we have ||m'(t)|| <t, as well as ||TO(i)|| <i^/2. 
The difference between jump times is denoted as 

The proofs of the asymptotic variance reduction make use of the following 
asymptotic expansions of the jump time differences of both processes : 
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Lemma 4.4. The difference between two jump times of the control process satisfies 



(4.3) 



and 



Lemma 4.5. The jump time variations of the process with internal state can he 
written in the following form (nl(e) being defined by ( 2.11\ ) : 

Ar„+i - AT,^+i +eAy„Vi + {0l+i+0n+i)O (e' + enl(e)) , (4.4) 

where 

6^ 



at" 



n+l 



^0 



m'(AT,^+i)ZT„ 
+ 0n+iO(e') 





n+l , a rniJ 



(4.5) 



and, correspondingly, 



An -6^6 ^^-+1'^ Zt 



-&^to(ATO+i)V^(XtJV„ 



1 



-b m 



'n+l 

"a7 



]VS{XTjVn+9n+,0{e'). 



(4.6) 



It is also useful to recall that according to ( 2.2 )-( 2.13 1, the following hold : 

|AT„+i| <C6'„+i, |Ar,j_|_;^| <C6'„+i. 

Finally, we will also need a different expansion of the jump times of the process 
with internal state : 



Lemma 4.6. When using nl2(e) as defined by (2.111, the jump time variations of 
the process with internal state can be written in the following form : 



Ar„+i = 



9n+l b'^T, 



^0 ^0 ^ 

+ ((?3^i+0„+i)O(enl2(e) + 



(Zt,.+i - ZrJ+eATob^nVSiXrjVn 



(4.7) 



Proof. The proof is based on Duhamel integration of the ODE (2.4) on [r„,r„+i], 
as in the proof of Lemma 4.7 of the companion paper Following equation (4.24) 
iniH, we get forte [r„,T„+i] : 



Zt = e-(*-^" " ^T„ + e (id -e- ' ' ) r, V5(XtJ V„ 



9l^,+e^+i)0{cse^ + cpe'+'). 



(4.8) 



Multiplying both sides by r^, recalling that Assumption 2.4 implies \\t^\\ = 0{e ^) 
and using Assumption 2.3 we get for t£ [T„,T„+i] : 



T,Zt=ree-(*-^")"'' 'ZT„+e(ld- 



-(t-T„)r.- 



(4.9) 
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Moreover, integrating (4.8 1 on [T„,Tn+i] yields 

+ eiT„+i-Tn)T,VS{XTjVn + i0l+i + el+^)O{cse^ + CFe'+^) 



(4.10) 



Evaluating equation (4.91 at t = Tn+i, add adding equation (4.101, we get 
Ztdt + T,ZT„^, =t,Zt„ +e(T„+i -r„)r,V5(XTjV„ 

The estimates in Lemma 4.8 from ^6 (i.e., \ATn+i~ AT,"^_^_J^\={ef^_^_J^+0n+l)O{e)) 
yields 

Ztdt + T,ZT,^^, ^t,Zt„ +eAr°+ir,V5(XTjV„ 



+ {el+i + en+i) O (cs(ei+* + 6^) + (c^ + l)ei+* + c^^ei+^)) . 
Finally, plugging the last equation in the estimate (|2.3|) yields the result. □ 



4.2. Analysis of variance of the coupling 

The present section is devoted to the proof of the following theorem: 

Theorem 4.7. Assume the assumptions of Section\2.^ hold, and that kd>l where 



k>2 is defined in (2.3). Then the difference between the process with internal state 

:0(e + e*'+nl2(e)), p>l, (4.11) 



and the coupling process satisfies 



where nl2(e) is defined in (2.11 1. 



The proof of Theorem 4.7 relies on a number of steps that will be detailed by a 
series of lemmas. We will make use of the following random sequences, defined as 
the position of both processes after n jumps, 

S„:=XT„=Xo + e^Ar^+iV^, n>0, (4.12) 



as well as 



= X|,e=Xo + e^Ar„ViV™, n>0. 



(4.13) 



We will also make use of the following random integers : 

Definition 4.8. The random integers iV> 1 and N'^> 1 are uniquely defined by: 
TN<t/e^<TN+i, T^c<i/e^<T^.+^. (4.14) 
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The first lemma bounds the difference at time t/e^ between the two processes 
Xf^^2 and X^^^2 by expressing it in terms of differences of positions and jump times 
of both processes after the same random number of jumps, S„ and 'E.'^ : 



Lemma 4.9. The difference between the rescaled process with internal state :- 
Xi/^2 and the rescaled coupling process X^"'^ :=X^^^2 satisfies: 



\XI-X'^''\<\EN~Ej,\+e\TN-n\ + i0N+i+ON^+i)Oie) 

< I Satc ~ S$vc I + e I T^c -T^.\ + {9n+i + ^tv.+i ) 0(e) 



(4.15) 
(4.16) 



Only one of the two estimates (4.15 1-(4.16 1 is necessary in the remainder of the 
proof, but we detail both to highlight the symmetry. 



Proof. By definition, we have 



XI = EN + e{t/e^-TN)VN, 
X'f'=E'},.+eii/e^-T^.)VN, 



so that by Definition 4.8 of N and N"^, and by realizing that \t/e —T^a 
and |i/e^ — TtvI <C0Ar+i, we get 



<ce 



X} 



'^i I 



-ATI 



\xt-xr\ 



< s 



(a) 

I Sato — S7v| + (6'Ar+l +0Ar^ + l)O(e). 



To analyze (a), we consider three cases: 

N^Nc. Then (a)=0. 
N>Nr. Then 



(a) 
e 



(4.17) 



N-l 

71 = N a 

^ rpc _rpC _ ATiC 4_T'^ — T*^ 

<T^ — T]\[ + C6n<:+i, 

where in the last line we have used that, AT^c_^_i < COn^+i^ and, by defi- 
nition, Tn <t/e'^ kT^c^i- 
Nr>N. Then 



(a) 
e 



E 



-T] 



N 



<7'7v+i — Tat — AT/v+i +Tjv — < Tn —T^ + C9 



N+l, 



where in the last line we have used that by definition T^a <t/e^ ^TV+i. 
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This proves (4.151. By symmetry, we get: 
(b) 



which yields (4.16). 



□ 



In the next lemma, we estimate the supremum of the difference of the position 
of both processes after n jumps for n£ [0,nc]. 

Lemma 4.10. Let > 1 be a deterministic integer verifying — 0{e^^). Then 

i/p 

= 0(e + e*+nl(e)). (4.18) 



e( sup 



Proof. First, using the estimates of jump times (4.3)-(4.4), we can decompose the 
differences of positions of both processes as follows : 

n-l 

m— 
n-l 

= ttm+l +'^m+l +''m+l, 



where, by definition, 

f^m+l '■— e^"^ (1 ~ ^m+l) ^e(>^m)"^"'^m 

+£2^ (Ao6^m (^) V5(S^) - 
+e6^m'(AT0+i)ZT,„V„„ 
l^r„+i:-((?^+i+ft„+i)0(e3 + eMe))- 
Since is Lipschitz, we have that am+i < C'e^ |Sm — and we can write 



for n < , which yields 



|S„-5^|< sup 
yo<i<ne 

Note that 



m=0 



- sup 

0<i<nj 



m— 1 



E 



£|r,„| (l + Ce2 + ... + (Ce2D,. (4.19) 



(l + Ce2 + --- + (Ce2)"')<C. 

Now, we can remark that the discrete time process {AIi)^^^ := (^X]L=i^m) is a 
martingale with respect to the filtration Ti = (7{{9m,Vm~i)\l<m<l) with A/o = 0, 
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SO that we can apply the Burkholder-Davies-Gundy upper bound (which is a simple 
consequence of Doob maximal inequality here, see^, for some p>l: 



e( sup \MA<CpnP/^-'J2^m~Mi_,n. 



In the present case, since \Zt^ \ — 0{e^): 
Moreover, 



El 



El 

m— 1 



m— 1 



<C(e + nl(e)f , 



and the result follows from (4.19). 



(4.20) 



□ 



In the same way, we can estimate the difference of the n^-th jump time of both 
processes : 



Lemma 4.11. Let > 1 a deterministic integer. Then 

efif sup \Tn-T:i\A <Cp{e' + nhie)) 

\0<n<nt J 



(4.21) 



Proof. First, using the estimates of jump times (4.3l-(4.7), we can decompose the 
differences of jump times of both processes as follows : 

ri-l 
n-1 



m=0 



where, by definition. 



(0i+i%„+i)eO(e'+nl2(e)). 



L 'm+l 



We can write for n<ne: 



\Tn~K\< 



Ao 



- sup 

0<i<n, 



-El 

m— 1 
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Now, we can remark that the discrete time process (M;)^^^ :— ^X)rn=i'^"') ^ 
martingale with respect to the filtration J-j = tr ((6'm, Vm_i)|l < rn< Z) with Mo = 0. 



Using the Burkholder-Davies- Gundy inequahty ( |4.20p , together with: 
yields 



E(|m™^)<CeP^ 



E sup 

\ 0<l<n, 



<Ce 



5-1 



Using: 



An 



{Zt^~Zt„) 



and plugging in the rest term, we finally get (4.21). 



□ 



We can now conclude with the proof of Theorem 4.7 



Proof. [Proof of Theorem 4.7 We consider a given integer rt^ > 1. We decompose 
the difference between both processes as 



(4.22) 



and use (4.15) to write 



\XI-Xl''\\N<n,< sup |S„-S^|+e sup \T^-T^\ + {eN+i+eN^+i)0{e). 



0<n<n, 



Using now the estimates (4.18) and ( |4.21[ ) in (4.23), we get 
E ( I X| - 1^<„ J < C(e + + nh (e)) 

= 0{e + e' + nh{e)), 



(4.23) 
(4.24) 
(4.25) 



where in the last line we have used a technical lemma (Lemma Appendix A.l ). 

It remains to control the probability of the event {N>n^}. To this purpose, we 
consider i i— >■ TVj S N the Poisson process uniquely defined by: 

Nt Nt + 1 

n—1 n—1 

By definition of N and Nc, there is a constant C such that N(ji/^2 > N and A^ct/e^ ^ 
A^'^, so that we can write: 



|A|-A?''| <eC 




Ct/e2>ne 



(4.26) 
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Let US choose n^>4:eCt/e^, a choice that satisfies ne=C'(e^^). A second technical 



lemma (Lemma Appendix A. 2) impUes that 



l-N>n. 



)<eCp2-i/(^^'), 



(4.27) 



so that the contribution of the latter to the final estimate (4.11 ) is negligible, since 
it is exponentially small with respect to e~^. The proof is complete. □ 



4.3. Discussion about the sharpness of the coupling 

One may ask about the "sharpness" of the precise estimates of the coupling, given 



mainly by Theorem |4.7[ (see also (4.18) (position shift), and (4.21) (time shift)). 
Let us discuss the case where cp = 0, Sg [l/fc,l/(fc — 1)]. In this case the dominant 
error term in the jump times expansion of the internal state model in Lemma 4.8 
in 1211 is due to the non-linearity of the turning rate: cxe'^^~^- Since this term is due 
to the internal state mechanism, it can be conjectured that the difference between 
the density of the internal state model and gradient sensing model will be at least 
of this order. On the other hand, the coupling estimate in Theorem |4.7| is controlled 
by the same term: cxe''^~^ , suggesting the sharpness of the latter. 



5. Numerical illustration 



In this section, we demonstrate the validity of the analysis above. For the numerical 
experiments, we restrict ourselves to one space dimension. For the process with inter- 



nal state, we use internal dynamics that are given by the scalar cartoon model ( 2.5 ). 
In all experiments, r is chosen independently of e so that 6=1. The corresponding 



turning rate is given by (2.6) with Ao = l. For the control process, we choose the 



linear turning rate (2.14), with parameters (2.22) such that the two processes have 



the same diffusion limit; in particular, 6 = 1. The physical domain is x€ [0,20], and 
we use reflecting boundary conditions, i.e. the bacterial velocity is reversed when 
a: = or a; = 20. We fix a scalar bimodal chemoattractant concentration field 



S{x) = a (exp (x - C)') + cxp (-/3 (x -??)')) , 
in which a = 5, /? = 1, C = 7.5 and rj = 12.5. 



(5.1) 



Single bacterium as a function of time. In a first experiment, we simulate 



a single bacterium evolving according to the fine-scale model (2.4), in which we 
set the parameters to e = 0.2 and r = l. We compare this evolution to that of a 



bacterium that satisfies the corresponding control process (2.12). Both simulations 



are performed using the same random numbers starting from the initial condition 
Xq^Xq=8 and Vb = Fp'^ = -1-1; the bacterium with internal state starts from = 
S{Xo). The time step 6t = 0.1. The results are shown in figure [I We see a very good 
coupling initially, which degrades over time. Note the time shift in the short time 
picture. In the long time picture, coupling is completely lost at some point. 



November 23, 2011 6:7 WSPC/INSTRUCTION FILE 



variance 



24 



Variance reduced simulation of chemotaxis 
14 




Fig. 2. Evolution of a single bacterium evolving according to the fine-scale process ( [2^ (solid line) 
and the corresponding control process | |2.12| (dashed) on short (left) and long (right) time-scales. 




^ 10" 



10" 



10-* 



T = IO.O/e^ 
T = 20.0/£^ 
T = 30.0/e^ 
theory 



0.01 



0.02 



0.05 



0.1 



Fig. 3. Empirical mean (left) and variance (right) of the difference between the fine-scale process 
| |2.4| and the corresponding control process ( |2.12[ l as a function of e for different values of the 
reporting time T. The theoretical slope is indicated with a dashdotted line. The sample size is 
AT = 10000. 



Expectation and variance as a function of e and t. Next, we repeat the 
experiment using N = 10000 particles and compute the empirical mean and variance 
of the coupling, i.e., 

= resp., = . 

i i 

As fine-scale parameters, we choose t = 1, Xo = b=l and several values of e. The 



chemoattractant concentration is again given as (5.1), now with a~f3 — l, ^ = 7.5 
and 7^ = 12.5. The time interval for the computation is tG [0,30/e^]. Figure [s] shows 
the dependence in e of the coupling, by plotting the empirical mean and variance 
defined above as a function of e for different values of the reporting time T. The 
results shown in figure [3] are in clear accordance with the theoretical slope predicted 
by the asymptotic analysis. Figure [4] shows the mean and variance of the coupling 
difference as a function of time. The time dependence of the variance has not been 
analyzed mathematically in Section [4j the specific behaviour viewed in figure |4] and 
is probably due to: (i) sufficiently short diffusive times; (ii) the specific double-well 
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>5, 




Fig. 4. Evolution of th e empirical mean (left) and variance (ri ght) o f the difference between the 
fine-scale process | |2.4| l and the corresponding control process as a function of t = t/e^ for 

different values of e. The sample size is Af = 10000. 





T = 10.0/e^ 
T = 20.0/e^ 
T = 30.0/e^ 



Fig. 5. Empirical mean (left) and variance (right) of the difference between the fine-scale process 
| |2.4[ l and the corresponding control process l |2.12| as a function of r for different values of the 
reporting time T. The sample size is A'^ = 10000. 



form of the chcmoattractant potential. 

Expectation and variance as a function of r. Finally, in a last experiment, we 
illustrate the dependence on t. To this end, we again simulate A'^= 10000 particles 
choosing Aq = 1, e = 0.1 and Xq = 7.5, for different values of r. The results in figure[5] 
show that the variance quickly increases with r until it reaches a plateau for t > 1. 



6. Simulation on diffusive time scales 



In this section, we consider a simulation of the density of an ensemble of particles, 
with and without variance reduction and/or reinitialization, as described in Sec- 
tion 3.1 We again restrict ourselves to one space dimension, with domain xG [0,20] 



and periodic boundary conditions. In this case, the kinetic equation corresponding 
to the control process reduces to the system 



dtp". - 



dxP+ = - 



2 

X^ix-l) 



-P- 



(6.1) 



-P- 
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Fig. 6. Bacterial density as a function of space at f = 50/e^ witiiout variance reduction. Left: one 
realization. Rigiit: mean over 100 realizations and 95% confidence interval. The solid line is the 
estimated density from a particle simulation using the process with internal state; the dashed line 
is estimated from a particle simulation using the control process. Both used A'^ = 5000 particles. 
The dotted line is the solution of the deterministic PDE ( |2.20| . 



of two PDEs, which is straightforward to simulate using finite differences. 



We fix the chemoattractant concentration field as (5.1 1, with parameters a 



(3 = 1, ^ = 7.5 and 77 = 12.5. For the internal dynamics, the same model ((2.5)-(2.6l 



(2.14|)-( |2.22| )) is used. The parameters are e = 0.5, Ao = l, t = 1, St = 0.1. 

All simulations are performed with A^ = 5000 particles. The initial positions are 
uniformly distributed in the interval a; €[13,15]; the initial velocities are chosen 
uniformly, i.e., each particle has an equal probability of having an initial velocity of 
±e. The initial condition for the internal variable is chosen to be in local equilibrium, 
i.e., Yq ^ S{Xq). The initial positions and velocities of the control particles are 
chosen to be identical. 



We discretize the continuum description (6.1) on a mesh with Ax = 0.1 using a 
third-order upwind-biased scheme, and perform time integration using the standard 
fourth order Runge-Kutta method with time step Stpde ~ 10^^. The initial condition 
is given as 

+ ( n^ -( 0^ [13,15], 
p+(a;,0)=p {x,{)) = \ . (6.2) 

10, otherwise. 

Simulation without variance reduction. First, we simulate both stochastic 
processes up to time i = 50 (t = 50/e^) and estimate the density of each of these 
processes nN{x,t), resp. n^(a:,f), without variance reduction. The density is ob- 
tained via binning in a histogram, in which the grid points of the deterministic 
simulation are the centers of the bins. Figure [6] (left) shows the results for a sin- 
gle realization. We see that, given the fluctuations on the obtained density, it is 
impossible to conclude on differences between the two models. This observation is 
confirmed by computing the average density of both processes over 100 realizations. 
The mean densities are shown in figure |6] (right), which also reveals that the mean 
density of the control process is within the 95% confidence interval of the process 
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Fig. 7. Bacterial density as a function of space at t = 50/e^ with variance reduction and reinitial- 
ization. Left: variance reduced density estimation of one realization with Af = 5000 particles (solid) 
and deterministic solution for the control process | [Z20l l (dashed). Right: mean over 100 realization 
and 95% confidence interval (solid) and the deterministic solution for the control process l |2.20| 
(dashed). 

with internal state. Both figures also show the density that is computed using the 
continuum description, which coincides with the mean of the density of the control 
particles. 



Simulation with variance reduction. Next, we compare the variance reduced 



estimation (3.1) with the density of the control PDE. We reinitialize the control 
particles after each coarse-scale step, i.e., each k steps of the particle scheme, where 
k6t = Stpde, (here k = l). The results are shown in figure [7] We see that, using this 
reinitialization, the difference between the behaviour of the two processes is visually 
clear from one realization (left figure). Also, the resulting variance is such that 
the density of the control PDE is no longer within the 95% confidence interval 
of the variance reduced density estimation (right figure). We see that there is a 
significant difference between both models: the density corresponding to the control 
process is more peaked, indicating that bacteria that follow the control process are 
more sensitive to sudden changes in chemoattractant gradient. This difference can 
be interpreted from the fact that the bacteria with internal state do not adjust 
themselves instantaneously to their environment, but instead with a time constant 



Simulation with variance reduction but without reinitialization. We also 



compare the variance reduced estimation (3.1 1 with the density of the control PDE 
without performing any reinitialization of the control process to restore the coupling. 
Figure [s] (right) shows the mean estimation of the density over 200 realizations, as 
well as the density of the control process. As evidenced by the error bars, simulating 
a single realization of the process, even with variance reduction, is not able to 
reliably reveal this difference. This phenomenon is due to the degeneracy of the 
coupling on long diffusive times. This is also illustrated in figure [s] (left), which 
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Fig. 8. Bacterial density as a function of space at t = 50/e^ with variance reduction. Left: variance 
reduced density estimation of one realization with A'^ = 5000 particles (solid) and deterministic so- 
lution for the control process ( |2.20| l (dashed). Right: mean over 200 realization and 95% confidence 
interval (solid) and the deterministic solution for the control process l|2.20|l (dashed). 
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Fig. 9. Bacterial density as a function of space at t = 50/e2 with variance reduction and reini- 
tialization. Left: variance reduced density estimation of one realization with Af = 5000 particles 
(solid) and deterministic solution for the control process | |2.20| (dashed). Chemoattractant given 
by l |5.1| with = 5. Right: mean over 100 realization and 95% confidence interval (solid) and the 
deterministic solution for the control process <|2.20|l (dashed). 



compares the variance reduced density estimation of a single realization with the 
density of the control PDE. Note that, as predicted by the analysis, the variance is 
much larger in regions where VS'(a;) is large. 

Finally, we repeat the experiment with a larger chemoattractant gradient. We 
again choose the chemoattractant field as (5.1 ), with the same parameters as above, 
except that we now take a — 5. We also use the same discretization parameters as 
above. The variance reduced estimation (3.1 1, as well as the density of the control 
PDE, are shown in figure[9] We see that, although the difference between the process 
with internal state and the control process becomes much larger, the variance is still 
very well controlled by the coupling. 
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7. Conclusions and discussion 

In this paper, we studied the simulation of stochastic individual-based models for 
chemotaxis of bacteria with internal state. We have used a coupling with a simpler, 
direct gradient sensing model with the same diffusion limit to obtain an "asymp- 
totic" variance reduction. Throughout this work, we have assumed that the com- 
putation of the bacterial density using the kinetic description of the direct gradient 
sensing model can be performed accurately with a grid-based method. Note that in 
higher spatial dimensions, this may become cumbersome, and the grid-based com- 
putation should be carried out at the level of the advection-diffusion description 
(or any desired moment system) only, using a second level of coupling between the 
velocity-jump process associated with the kinetic description, and the stochastic 
differential equation associated with the advection-diffusion limit. This is left for 
future work. 

The current algorithm allows to explore the differences between the fine-scale 
process with internal state and the simpler, coarse process using a fixed number of 
particles, independently of the small parameter e. However, the computational cost 
to simulated an individual particle over diffusive time-scales becomes very large 
for e— s-0. In future work, we will therefore also study the development of truly 
"asymptotic preserving" schemes in the diffusion asymptotics, in the sense that 
the computational cost of the simulation of processes is independent of the small 
parameter e. This will require to deal with two kinds of difficulties: (i) We will need 
to use an asymptotic preserving method to solve the density evolution of the control 
model on the grid; and (ii) We will need to extrapolate forward in time the state of 
the fine-scale simulation. 

The first difficulty implies that, instead of solving the full kinetic equation as- 
sociated with the control process, we only solve its diffusion limit. This may imply 
the use of a second level of variance reduction, coupling the control velocity jump 
process and its limiting drift-diffusion process. 

The second difficulty, extrapolation in time, is related to the equation-free 1121221 
and HMM HEI types of methodologies; ideas of this type can be traced back to 
Erhenfest One approach is to use a coarse projective integration method EEH 
One then extrapolates the bacterial density over a projective time step on diffusive 
time scales, after which the projected density needs to be lifted to an ensemble of 
individual bacteria. For such methods, besides the effect of extrapolation on the 
variance of the obtained results, also the effects of reconstructing the velocities and 
internal variables have to be systematically studied. 
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Appendix A. Technical lemmas, used in proof of Theorem |4.7| 

Lemma Appendix A.l (Technical lemma). Let {dn)n>i a sequence of indepen- 
dent exponential random variables with mean 1, and (T'„)„>i a strictly increasing 
sequence of random times with Tq =0 such that: 

• For any m>l, the sequence {0n)n>m+i 'is independent of the past 

• There is a constant C such that for all n > 0, 

-^dn+l<Tn+l—Tn<C9n-\-l- (A.l) 

Let Nt a random integer such that: 

Then for any integer fc > 0, there is a constant Ck independent of t .such that: 

Proof. The result may seem rather intuitive, but the proof is slightly tricky. First, 
we consider sub-intervals of [0,t] of the form: 

for some P>1 and p = 1 , . . . , P. 
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Step (i). The first step consists in proving the fohowing estimate: 

P t 



p=0 



CP 



sup VP(r„ eJp 



(A.2) 



Indeed, 



+CSO 



Jp= 

Using the independence of 0„+i with r„, and the assumption (A.l), it yields, 



P-p 



and finally changing the sum index p to P — p-{-l yields (A.2). 

Step (ii). The second step consists in the following decomposition. Let I — [a,b] 
be a finite interval of M+, T_i — — oo. Then, 



+ 00 



+ 00 



+ CX3 



^P(T„ ei) = Y,[ ^{Tn e /,r„_i ^ /) ^ P(r„+i, . . . ,r„+„, e /|r„ e /,t„_i ^ /) j . 

(A.3) 



n=0 



n=0 



m=0 



The key is to write 

+00 +00 +00 

^ P(r„ e /) - ^ p(r„_ i ^ /, r„ e /) + ^ P(r„_ i , r„ e /) 



n=0 



n=0 



+00 



n=0 



+ C30 



ri=0 



n=0 



(ai) 



We wish to iterate the decomposition of (ao) to (ai), (02), etc... For this purpose, 
remark that 



+00 



) = ^p(T„,...,r„+„G/) 

+00 

= HTn e /)P(T„+i , . . . , T„+,n e I\Tn G /) 



/ +00 



< ( ^P(T„G/)jP(0i + --- + 0™<C(6~a)) 
>0, 
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SO that we get in the end 



+00 +00 



(ao) = 5m ^^^^-^ iI,Tn,. . .r„+„ € /). 



Factorizing P(r„ e /,T„_i ^ /) with Bayes' formula yields (A. 3 1. 
Step (iii). The estimate (|A.2[) yields 



+00 



+ 00 



'1 



< 



t 

CP 



(A.4) 



p— m— 

Indeed for 1= [a, 6], using the independence of (^n+i, • • • ,^n+m) from (Tn,Tn-i): 

+ 00 +CX3 

^P(r„+i,...,r„H.„e/|r„e/,T„_i^/)<^p(0i + ---+^™<C(6-a)), 



m=0 



m=0 



and 



+00 



^ P(T„ e /,r„_i ^ /) = P(37i > 0|T„ e /) < 1. 



The two last estimates in (A.3) yields: 
+00 



+■00 



sup VP(T„e/p)<Vp 0i + --- + 0,„< 

^ n— m— 



and (A.4) follows. Finally, we choose P= [t + lj so that the right hand side of (A.4 1 
does not depend on t any more. Finally we use the exponential convergence of the 
terms of the two series in the right hand side of ( A.4[ ) (using for instnace large 
deviation estimates for the second) one, to conclude that the two series indeed 
converge, and conclude the proof. □ 



Lemma Appendix A. 2 (Technical lemma). Let {On)n>o °' sequence of ex- 
ponential random variables of mean 1, and ti-^ NtGN the Poisson process uniquely 
defined by: 

Nt Nt + 1 

^0n<t< J2 
n— 1 n — 1 

Let uqCzN such that nQ>Aet, pi>l, and p2>^- Then there is a constant Cp-^^p^ 
independent of (riQ^t) such that: 



E 



Nt + l 
n=l 



P2 



LA't>no 



Proof. 
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Step (i). 

Let us condition on the different possible values of Nt>no: 



E 



Nt + l 
n=l 



P2 



+00 



LjVt>no 



m— riQ + l 

+00 

E E| 



E^ 
E^ 



P2 



P2 



Nt=m ]V{Nt=m) 



Nt=m e 



m! 



Step (ii). 

Let us denote (Pt,m,]Et,m) = (P( ■ =m-) ,E(- ='m)) the probability /expectation 
conditionally on the event {Nt=m]. Conditionally on {Nt = ra], the se- 
quence (01, . . . ,6'„j+i) has the same distribution as (tC/(i) ,i[/(2) j ■ ■ ■ 5it^(m) ~ 
tJ7(m_i),i — +6'), where (t7(i) <•••< is the order statistics of m inde- 
pendent random variables (Ui, . . . ,Um) uniformly distributed on the interval [0,1], 
and 9 is an independent exponential random variable of mean 1 (see Section 2.4 
of 1221 for classical properties of Poisson processes). Thus denoting C/(o) =0 we can 
write: 



E 



Nt + l 

E ^^'"^ 

m— no + l 



^Nt>no 



,P1 



n=l 



(A.5) 



Using the inequality {a + b)^ <CpaP + for any a,6>0, Jensen inequality, and de- 
noting J7(m+i) — 1, there is a constant Cp-^^^p^ such that: 



pi 



< 



/ m+l 

Cp,,p2 0''^fi+iPif^(m-f ir-1 ^ ([/(„)- 



\PlP2 



(A.6) 



Step (iii). 

Let us compute Et_m (([/(„) — U(^n-i)Y^^^) for ?i = 1, . . .,m+l. First, let us recall the 
probability density of a couple [U ({-^ ,U (^j)) for 1 <i < j <m in an order statistics of 
size to: 

P([/(i) e [uj,Ui-|-dui],C/(j) e [uj,Uj + duj]) 

^m\V{Ui<--- <Ui-i<Ui)xP{Uie[ui,u., + dui\) 

X P(uj < t/j+i < • • • < Uj-i < Uj) X P(C/j e +dMj]) X P(mj < < • • • < Urn) , 
(u,-ii, )■'"'"' (1-u,)"-^' 



- to! ^ — - — TT ^ — — — — dui dui . 

{m-j)\ ' 



(A.7) 
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As a consequence, 

Et.™((f/(„)-C/(„-i))'"''') = 



0<M„_i<n„<l 

m!(pip2)! 



■(n-2) 



{Un-Un-lY 



(m — n)! 



-dun^i dun 



0<u„_i<u„<l 



(n-2)! (P1P2)! (m-n)! 



'7i!(piP2)! 

(m+pip2)!' 



where in the last hne we have used (A. 7) with appropriate {i,j,m) to compute the 
integral. In the same way, the probability density of C/(,j) for l<i<m in an order 
statistics of size m is given by: 

fl — w-V"^* 



so that 



Et.^(c/fif^)=m!(pip2)! 

_ m!(pip2)! 



(i-1)! (m-i)! 

^'^^ (l-m)"'-i 
o<ui<i (P1P2)! (m-1)! 



(rn+PiP2)!' 



and in the same way. 



Et,™((l-C/(„,)r 



"T^!(PlP2)! 

(m+pip2)!' 



Step (iv). 

Remarking that Em(^'^^^^) = (P1P2)!, and inserting the computations of Step (ii) 
in (|X5l)-(rO) yields: 



E 



Nt + l 
n=l 



P2 



< 



+00 



m— 710 + 1 

Using the inequalities ^"l^+plp^y' < 1, e"* < 1 and 

-^-^ m! ~ rin! ' 

m— no + l 



it yields 



E 



Nt + l 
n=l 



P2 



\ frio+piP2 
'-JVt>no j ^CpiP2 ^^^^^1 
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SO that in the end, using the assumption no > 4et and the Stirling formula, 

P2 \ 



El 



Nt + l 

n=l 
<Cr, 



^0 



which yields the result. 
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